home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / mathpack / daxpy.f < prev    next >
Text File  |  1989-08-17  |  4KB  |  99 lines
  1. c   imsl routine name   - vbla=daxpy                                    vbdb0010
  2. c
  3. c-----------------------------------------------------------------------
  4. c
  5. c   computer            - vax/double
  6. c
  7. c   latest revision     - january 1, 1978
  8. c
  9. c   purpose             - compute a constant times a vector plus
  10. c                           a vector, all double precision
  11. c
  12. c   usage               - call daxpy (n,da,dx,incx,dy,incy)
  13. c
  14. c   arguments    n      - length of vectors x and y. (input)
  15. c                da     - double precision scalar. (input)
  16. c                dx     - double precision vector of length
  17. c                           max(n*iabs(incx),1). (input)
  18. c                incx   - displacement between elements of dx. (input)
  19. c                           x(i) is defined to be..
  20. c                           dx(1+(i-1)*incx) if incx.ge.0 or
  21. c                           dx(1+(i-n)*incx) if incx.lt.0.
  22. c                dy     - double precision vector of length
  23. c                           max(n*iabs(incy),1). (input/output)
  24. c                           daxpy replaces y(i) with da*x(i)+y(i) for
  25. c                           i=1,...,n.
  26. c                           x(i) and y(i) refer to specific elements
  27. c                           of dx and dy, respectively. see incx and
  28. c                           incy argument descriptions.
  29. c                incy   - displacement between elements of dy. (input)
  30. c                           y(i) is defined to be..
  31. c                           dy(1+(i-1)*incy) if incy.ge.0 or
  32. c                           dy(1+(i-n)*incy) if incy.lt.0.
  33. c
  34. c   precision/hardware  - double/all
  35. c
  36. c   reqd. imsl routines - none required
  37. c
  38. c   notation            - information on special notation and
  39. c                           conventions is available in the manual
  40. c                           introduction or through imsl routine uhelp
  41. c
  42. c   copyright           - 1978 by imsl, inc. all rights reserved.
  43. c
  44. c   warranty            - imsl warrants only that imsl testing has been
  45. c                           applied to this code. no other warranty,
  46. c                           expressed or implied, is applicable.
  47. c
  48. c-----------------------------------------------------------------------
  49. c
  50.       subroutine daxpy  (n,da,dx,incx,dy,incy)
  51. c
  52. c                                  specifications for arguments
  53.       double precision   dx(1),dy(1),da
  54.       integer            n,incx,incy
  55. c                                  specifications for local variables
  56.       integer            i,iy,m,mp1,ns,ix
  57. c                                  first executable statement
  58.       if (n.le.0.or.da.eq.0.d0) return
  59.       if (incx.eq.incy) if (incx-1) 5,15,35
  60.     5 continue
  61. c                                  code for nonequal or nonpositive
  62. c                                    increments.
  63.       ix = 1
  64.       iy = 1
  65.       if (incx.lt.0) ix = (-n+1)*incx+1
  66.       if (incy.lt.0) iy = (-n+1)*incy+1
  67.       do 10 i=1,n
  68.          dy(iy) = dy(iy)+da*dx(ix)
  69.          ix = ix+incx
  70.          iy = iy+incy
  71.    10 continue
  72.       return
  73. c                                  code for both increments equal to 1
  74. c                                    clean-up loop so remaining vector
  75. c                                    length is a multiple of 4.
  76.    15 m = n-(n/4)*4
  77.       if (m.eq.0) go to 25
  78.       do 20 i=1,m
  79.          dy(i) = dy(i)+da*dx(i)
  80.    20 continue
  81.       if (n.lt.4) return
  82.    25 mp1 = m+1
  83.       do 30 i=mp1,n,4
  84.          dy(i) = dy(i)+da*dx(i)
  85.          dy(i+1) = dy(i+1)+da*dx(i+1)
  86.          dy(i+2) = dy(i+2)+da*dx(i+2)
  87.          dy(i+3) = dy(i+3)+da*dx(i+3)
  88.    30 continue
  89.       return
  90. c                                  code for equal, positive, nonunit
  91. c                                    increments.
  92.    35 continue
  93.       ns = n*incx
  94.       do 40 i=1,ns,incx
  95.          dy(i) = da*dx(i)+dy(i)
  96.    40 continue
  97.       return
  98.       end
  99.